capture program drop ccn_fit
program define ccn_fit, rclass
	local model_type = "`1'"
	local inv_var = "`2'"
	local dep_var = "`3'"
	local raw_control_list = "`4'"

	preserve
	local minlevel = 0.001
	local controls = "`raw_control_list'"
	if ("`model_type'" == "tobit") {
		ren `inv_var' I
		ren `dep_var' S
		qui sum I
		local samp = r(N)
		noisily qui tobit I, ll(0) vce(robust) noconstant
		qui predict hat
		qui matrix coeff_mat = e(b)		
		qui matrix coeff_Vmat = e(V)		
		local muq=coeff_mat[1,1]
		local sigma = coeff_Vmat[2,2]^(1/2)
		*X-quantiles
		qui pctile mu = I, nq(100) genp(q)
		*Getting muq_hat_star based on Tobit
		gen qq=q/100
		gen Z=invnormal(qq)
		gen muH_Tobit=`muq' + `sigma' * Z
		
		drop qq
		keep q mu muH_Tobit 
		order q mu muH_Tobit 
		qui keep if q != .
		save "$datapath/002_Fit.dta", replace //this file is saved to construct the plot of muhat_star across all methods
	}
	else if ("`model_type'"=="bilog") {
		ren `inv_var' I
		ren `dep_var' S
		gen cens_ind = (I==0)
		gen q=-_n if _n<=5000
		cumul I, generate(CDF_I)
		sum cens_ind
		local F0 = r(mean)
		local band = $band
		reg CDF_I I if I>0 & I<`band'
		local f0lim = _b[I]
		gen CDFmuH_bilog_lb=1 - (1-`F0')*exp(-(`f0lim'/(1-`F0'))*q)
		gen CDFmuH_bilog_ub=`F0'*exp((`f0lim'/`F0')*q)
		keep q CDFmuH_bilog_lb CDFmuH_bilog_ub
		qui keep if q != .
		append using "$datapath/002_Fit.dta"
		save "$datapath/002_Fit.dta", replace //this file is saved to construct the plot of muhat_star across all methods
	}

	else if ("`model_type'"=="nohighpeaks") {
		ren `inv_var' I
		ren `dep_var' S
		gen cens_ind = (I==0)
		gen q=_n if _n<=100
		gen qq=q/100
		cumul I, generate(CDF_I)
		sum cens_ind
		local F0 = r(mean)
		local band = $band
		reg CDF_I I if I>0 & I<`band'
		
		local slope = _b[I]
		gen muH_NoHighPeaks=(qq -`F0')/`slope'
		drop qq
		keep q muH_NoHighPeaks
		qui keep if q != .
		merge 1:1 q using "$datapath/002_Fit.dta"
		drop _merge
		save "$datapath/002_Fit.dta", replace //this file is saved to construct the plot of muhat_star across all methods
	}
	
	
	else if ("`model_type'" == "heterotobit") {
		ren `inv_var' I
		ren `dep_var' S
		qui sum I
		local samp = r(N)
		gen cens_ind = (I==0)
		*Run Tobit
		qui tobit I, ll(0)
		qui predict hat
		qui matrix coeff_mat = e(b)
		local sigma = coeff_mat[1,2]^(0.5)
		*X- quantiles
		qui pctile muq = hat, nq(100) genp(q)
		*Getting muq_hat_star based on HeteroTobit
		gen qq=q/100
		gen Z=invnormal(qq)
		gen muH_HeteroTobit=muq + `sigma' * Z
		drop qq
		keep q muH_HeteroTobit
		qui keep if q != .
		merge 1:1 q using "$datapath/002_Fit.dta"
		drop _merge
		save "$datapath/002_Fit.dta", replace //this file is saved to construct the plot of muhat_star across all methods
	}
	else if ("`model_type'" == "heterounif") {
		ren `inv_var' I
		ren `dep_var' S
		qui sum I
		local samp = r(N)
		gen cens_ind = (I==0)
		gen q=_n if _n<=100
		*Getting muq_hat_star based on HeteroUnif
		gen qq=q/100
		qui sum I if I>0
		local pos_mean = r(mean)
		sum cens_ind
		local bunching = r(mean)
		local slope = ((1-`bunching')/2)/`pos_mean' 
		gen muH_HeteroUnif=(qq -`bunching')/`slope'
		drop qq
		keep q muH_HeteroUnif
		qui keep if q != .
		merge 1:1 q using "$datapath/002_Fit.dta"
		drop _merge
		save "$datapath/002_Fit.dta", replace //this file is saved to construct the plot of muhat_star across all methods
	}
	else if ("`model_type'" == "heterosymm") {
		ren `inv_var' I
		ren `dep_var' S
		qui sum I
		local samp = r(N)
		gen cens_ind = (I==0)
		qui gen pvec = .
		qui sum cens_ind
		scalar P = r(mean)
		qui replace pvec = P
		qui replace pvec = 0 if pvec < `minlevel'
		qui pctile muq = I, nq(100) genp(q)
		*need to keep track of which are dropped
		local skip_ind = 0
		qui sum muq
		if r(max)==0 {
			local skip_ind = 1
			qui drop muq
			}	
		qui tempfile main
		qui save "`main'", replace
		
		qui keep if q != .
		qui keep muq q
		if (`skip_ind' == 0) {
			*case 1: median not censored (heterosymm)
			qui gen muq_symm2 = .
			if (muq[100/2]>0) {
				qui gen muq_refl = .
				sort q
				qui replace muq_refl = -muq[100-_n] + (2*muq[100/2])
				qui gen muH_HeteroSymm = muq_refl if muq<=0
				qui replace muH_HeteroSymm = muq if muq>0
				qui drop *refl
			}
			*case 2: regression style (heterotobit)
			if (muq[100/2] == 0) {
				gen muqFit=0
				gen Z=invnormal(qq)
				replace muqFit=muq + `sigma' * Z
				if _rc==0 {
					qui gen muH_HeteroSymm = muqFit if muq == 0
					qui replace muH_HeteroSymm = muq if muq> 0
					qui drop muqFit
				}
				if _rc != 0 {
					local skip_ind = 1
				}
			}
		}
		keep q muH_HeteroSymm
		qui keep if q != .							
		merge 1:1 q using "$datapath/002_Fit.dta"
		drop _merge
		save "$datapath/002_Fit.dta", replace //this file is saved to construct the plot of muhat_star across all methods
	}
	else if ("`model_type'" == "heterosymmetry") {
		ren `inv_var' I
		ren `dep_var' S
		qui sum I
		local samp = r(N)
		gen cens_ind = (I==0)
		qui gen pvec = .
		qui sum cens_ind
		scalar P = r(mean)
		qui replace pvec = P
		qui replace pvec = 0 if pvec < `minlevel'
		qui pctile muq = I, nq(100) genp(q)
		*need to keep track of which are dropped
		local skip_ind = 0
		qui sum muq
		if r(max)==0 {
			local skip_ind = 1
			qui drop muq
			}	
		qui tempfile main
		qui save "`main'", replace
		
		qui keep if q != .
		qui keep muq q
		if (`skip_ind' == 0) {
			*case 1: median not censored (heterosymm)
			qui gen muq_symm2 = .
			if (muq[100/2]>0) {
				qui gen muq_refl = .
				sort q
				qui replace muq_refl = -muq[100-_n] + (2*muq[100/2])
				qui gen muH_HeteroSymmetry = muq_refl if muq<=0
				qui replace muH_HeteroSymmetry = muq_refl if muq>0 & muq<muq[100/2]
				qui replace muH_HeteroSymmetry = muq if muq>0 &  muq>=muq[100/2]				
				qui drop *refl
			}
			*case 2: regression style (heterotobit)
			if (muq[100/2] == 0) {
				gen muqFit=0
				gen Z=invnormal(qq)
				replace muqFit=muq + `sigma' * Z
				if _rc==0 {
					qui gen muH_HeteroSymmetry = muqFit if muq == 0
					qui replace muH_HeteroSymmetry = muqFit if muq> 0
					qui drop muqFit
				}
				if _rc != 0 {
					local skip_ind = 1
				}
			}
		}
		keep q muH_HeteroSymmetry
		qui keep if q != .							
		merge 1:1 q using "$datapath/002_Fit.dta"
		drop _merge
		save "$datapath/002_Fit.dta", replace //this file is saved to construct the plot of muhat_star across all methods
	}
	restore	
end

